Dynamics of a quantum phase transition in a ferromagnetic Bose-Einstein condensate 
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We discuss dynamics of a slow quantum phase transition in a spin-1 Bose-Einstein condensate. 
We determine analytically the scaling properties of the system magnetization and verify them with 
numerical simulations in a one dimensional model. 



Studies of phase transitions have traditionally focused 
on equilibrium scalings of various properties near the crit- 
ical point. Dynamics of the phase transition presents new 
challenges and there is a strong motivation for analyzing 
it. Nonequilibrium phase transitions may play a role in 
the evolution of the early Universe Their analogues 
can be studied in the condensed matter experiments. The 
latter observation led to development of the theory based 
on the universality of critical behavior which in turn 
resulted in a series of beautiful experiments Q . The re- 
cent progress in the cold atom experiments allows for 
time dependent realizations of different models under- 
going a quantum phase transition (QPT) 0, 0]. These 
experimental developments are only a proverbial tip of 
the iceberg, but they call for an in-depth theoretical un- 
derstanding of the QPT dynamics. 

A QPT is a fundamental change in ground state (GS) 
of the system as a result of small variations of an exter- 
nal parameter, e.g., a magnetic field [6]. It takes place 
ideally at zero absolute temperature, which is in strik- 
ing contrast to thermodynamical phase transitions. The 
most complete description of the QPT dynamics has been 
obtained so far in spin models 0, that are exactly solv- 
able. In these systems the gap in the excitation spectrum 
goes to zero at the critical point, which precludes the adi- 
abatic evolution across the phase boundary. It leads to 
creation of excitations whose density and scaliiig with a 
quench rate follow from a quantum version 7, 9] of the 
Kibblc-Zurek (KZ) theory 1, 2]. 

We study dynamics of a ferromagnetic condensate of 
spin-1 particles [lO|. For simplicity, we consider ID ho- 
mogeneous (untrapped) clouds: atoms in a box as in the 
experiment ll| with spinless bosons. We adopt the pa- 
rameters for our ID model such that the length and time 
scales are comparable to experimental ones [12.] . Assum- 
ing that the system is placed in a magnetic field B aligned 
in the z direction, one gets the following dimensionless 
mean-field energy functional [l^ 



Em = dz 



1 d^-t d^ 



2 dz dz 

IB 



Co 
2 



(1) 



where ^p-^ ~ (-01, -00, V'-i) describes the m — 0,±1 con- 
densate components, / dz^^ IV'mP = 1, and Fx,y.z are 
spin-1 matrices [l^l- The first term in ([T]) is the kinetic 



energy, the second and the fourth term describe spin- 
independent and spin-dependent atom interactions re- 
spectively, the third term is a quadratic Zeeman shift 
coming from atom interactions with a magnetic field. 
For ^^Rb atoms considered here ci < 0, which results 
in an interesting phase diagram due to the competition 
between the last two terms in ([l}. Restricting analysis 
to zero longitudinal magnetization case, and introducing 

q^Q/{n\ci\), n = *t^ 

one finds a polar phase for ? > 2, described by ~ 
(0,1,0), and the broken-symmetry phase where 



for < (7 < 2. The freedom of choosing the x±i results 
in rotational symmetry of the transverse magnetization 
on the (x, y) plane. The transition between these phases 
can be driven by the change of the magnetic field B im- 
posed on the atom cloud, q ^ Q ^ B^ IJ], which was 
experimentally done in [5|. 

The dynamics of a QPT depends on the rate of quench 
driving the system across the phase boundary. For very 
fast "impulse" transition, the system has no time to ad- 
just to the changes of the Hamiltonian and arrives in a re- 
gion where a new phase is expected with the "old" wave- 
function untouched during the evolution. Slow transi- 
tions are different: the system has time to "probe" vari- 
ous broken symmetry "vacua" in the neighborhood of the 
critical point where it gets excited. We are interested in 
evolutions that are fast enough to produce macroscopic 
excitations of the system, but slow enough to reflect scal- 
ings of the critical region. By comparing analytical find- 
ings to numerical simulations for experimentally relevant 
parameters we provide the first complete description of 
QPT dynamics in a ferromagnetic condensate. 

Fast transitions were realized in the Berkeley experi- 
ment . The 3D numerical simulations closely following 
this experiment were reported in |14j . Analytical studies 
of the evolution after "impulse" quench were presented in 
lil lij . The paper of Lamacraft flBj also discusses dy- 
namics of non instantaneous transitions in 2D spinor con- 
densates focusing on analytical predictions on the growth 
of the transverse magnetization correlation functions. 

We start with a qualitative discussion. Consider- 
ing small perturbations around the GS of the broken- 
symmetry phase one finds three Bogolubov modes as in 



2 



|13l | where quantum fluctuations are studied. In the long 
wavelength limit (important for slow transitions) there 
is only one nonzero eigenvalue mode: the gapped mode 
having eigenenergy A ~ ■\/4 — (7^. Suppose now that we 
drive the system from polar to broken-symmetry phase. 
The system reaction time to Hamiltonian changes in the 
broken-symmetry phase is given by the inverse of the gap: 
r ~ -g- 0, 0. For example, when r is small enough the 
evolution becomes adiabatic so the system adjusts fast 
to parameter changes. Right after entering the broken- 
symmetry phase, the reaction time is large with respect 
to the transition time, A/^, and so the system under- 
goes the "impulse" evolution where its state is "frozen" . 
The gapped mode starts to be populated around the in- 
stant i after entering the broken- symmetry phase: the 
system leaves the "impulse" regime to catch up with in- 
stantaneous GS solution. This happens when the two 
time scales become comparable: 1/A{i) ^ A/^\^^^. We 
consider here transitions driven by 



q{t) = 2-t/TQ, 
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where tq is the quench time inversely proportional to the 
speed of driving the system through the phase transition. 
For slow transitions of interest here, tq 3> 1, we obtain 
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In the following we analyze dynamics induced by a 
linear decrease of q{t) The evolution starts from 

i < 0, i.e., in the polar phase, and ends at t = 2tq {q = 
0). Such q{t) dependence is achieved by ramping down 
the magnetic field as ~ — tjTq. The initial state 
is chosen as a slightly perturbed GS in the polar phase, 
^'^ ~ ((5'0i,l/\/I+(5-0o,'5'0-i), where I^Vml < 1/v^are 
random. We generate the real and imaginary part of b'^ra 
at different grid points with the probability distribution 
p[x) — exp(— x^/2(T^)/\/27rcr. We take a ~ 10^^ to start 
evolution closely to the polar phase GS. 

To find the full numerical solution within the mean- 
field approximation, we integrate three coupled nonlin- 
ear Schrodinger equations for the tpm condensates that 
can be easily obtained by the variation of ([T|). During 
evolution we look at the magnetization of the sample 

The transverse magnetization. A total transverse (to the 
magnetic field in the z direction) magnetization reads 



MT(t)= / dz[f^iz,t) + f^(z,t)] 



dz mT, 
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and is experimentally measurable. It disappears in the 
GS of the polar phase and equals (1 — g^/4)/i in the 
broken-symmetry GS. Its typical evolution is depicted in 
Fig. [T] We see there that nothing happens in the po- 
lar phase. The system starts nontrivial evolution in the 
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FIG. 1: (color) Main plot: numerical solution (black solid 
line) vs. static prediction (red dashed line). The arrow depicts 
direction of evolution. Inset (a): the same as in the main plot 
plus a numerically obtained solution of the linearized problem 
(green divergent line). Inset (b): numerical data vs. fit to 
TQ > 10 data only (see text for details). In the main plot 
TQ — 10 (see [1^ for units). 
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FIG. 2: The vectors represent {f^{z), fy{z)) x 10=*. Plot (a): 
snapshot at q{t = 2.81) — 1.72, i.e., at the first peak in AItL 
(see Fig. Plot (b): snapshot by the end of time evolution: 
q(t — 20) = 0. The results come from the same numerical 
simulation as in Fig. [T](see for units). 



broken-symmetry phase at a distance i/Tg after the criti- 
cal point was passed. The magnetization grows fast from 
that point until it exceeds the static prediction and starts 
oscillations with the amplitude decreasing in time. We 
consider slow transitions. Therefore, by the end of time 
evolution, when 9 = 0, the system is in the slightly per- 
turbed ferromagnetic GS: globally MtL « 1 (Fig. [T]) and 
locally L'^mT{z) w 1 (Fig. [3]). We can now ask: Does the 
scaling ([3]) hold? To find out we define arbitrarily t as the 
instant when MtL intersects 1%. A fit to numerics for 
TQ > 10 yields ln£= (0.056 ± 0.01) -t- (0.332 ±0.002) In tq 
which confirms prediction ([3]). This fit is presented in Fig. 

[TJi, where the gradual departure of the numerical data for 
1/3 

TQ < 10 from t ^ Tq indicates that tq 3> 1 or 37ms has 
to be taken for the observation of 1/3 exponent: quench 
has to be slow enough to reflect the critical dynamics. 
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FIG. 3: (color) Magnetization of the system a.t t — 2tq 
{q — 0). The dashed lines facilitate observation of extrema 
coincidences. Results come from the same simulation as in 
Figs, [nil see [ij] for units. 



In the GS configuration of the broken-symmetry phase 
the vector {fx, fy) can have arbitrary orientation, so in 
the dynamical problem considered here it is interesting to 
find out how is this rotational symmetry broken. When 
unstable evolution starts, spatial correlations in magne- 
tization appear (Fig. [2^). In the subsequent evolution 
these correlations evolve such that the correlation length 
increases: see Fig. [2b obtained by the end of time evolu- 
tion. This is a generic picture though the details depend 
on the quench time tq and initial state of the system. 
This behavior suggests creation of spin textures 13 ■ 
In our case, topological textures are spin configurations 
where the magnetization direction varies in space so that 
the kinetic energy term in |T]) is not minimized, but mag- 
netization magnitude follows closely a GS result. Such 
structures appear in ID when the first homotopy group 
of the vacuum manifold M is nontrivial, which happens 
here: Tri{M) — Z [l9|. These textures are characterized 
by the winding number, ^ / dzj^ATg{fj; + ify), which 
is not conserved. Indeed, it reads -1-1 in Fig. [2^, while 
by the end of that evolution (Fig. [2]d) it equals 0. 

Are different stages of this evolution experimentally 
observable? Let's look at tq = 10 case presented in Figs. 
[T]|3l The evolution from the phase boundary to the first 
peak in magnetization Mr (the q = point) takes 2.81 x 
37ms = 104ms {2tq = 740ms). Both these time scales 
are well within the reach of the experiment 

The longitudinal magnetization. Initially, fz{z) « so 
that / dzfz ~ 0. The conservation of the latter allows 
only for creation of a network of magnetic domains (non- 
topological structures with fixed fz sign) having opposite 
polarizations. The domains appear by the time when the 
system enters unstable evolution and the maxima of |/^| 
tend to move towards the minima of m^ (Fig. E]). More 
quantitatively, we performed Nr evolutions starting from 
different initial conditions, but fixed a. As in the exper- 
iment [si, we average over these runs to wash out shot- 
to-shot fiuctuations. In Fig. [4] we plot the mean domain 
size: ^ = J2i^z{i)/Nr, where i = l,...,Nr and ^z{i) is 
the mean domain size in the j-th run. As shown in Fig. 
12a, for t<iwe observe ^ ~ /(V^-q ^) as for Mrit). The 



domains are formed on a time scale oi^ t. A simple anal- 
ysis based on KZ theory [l|,0] suggests that their charac- 
teristic post-transition size, ^, should be roughly given by 

Jp^* dtvs {t) , where Vs {t) is a sound velocity. There are two 
sound modes in the broken-symmetry phase that prop- 
agate both spin and density fiuctuations 13]: the faster 
(slower) one has velocity ~ (~ •\/g|ci|). Putting any 
of these as Vg into the integral, and assuming tq ^ 1 for 
the slower mode, we get 
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This result correctly predicts the scaling property of 
the size of post-transition "defects" as is evident from 
the overlap of different curves in Fig. lU which shows 
up for Tq > 25 or 0.9s. Quantitatively, we define ^ 
as the value of ^ averaged over q £ [1/2,1] to wash 
out post-transition fiuctuations. A fit got us In^ = 
(-0.38 ± 0.03) -\- (0.30 ± 0.01) In tq, in good agreement 
with ([5]). The fit was done to tq > 30 data and is pre- 
sented in Fig. [lb which illustrates that smaller tq data 
gradually departs from 1 /3 scaling law. 

Now we focus on the analytical calculations provid- 
ing predictions about early stages of time-evolution. 
We assume that the wave-function stays close to 
the polar phase GS, = {S^pi{t),l/y/L + 

SiJjo{t), S^p-i{t)) exp(— ?/xt), where the chemical potential 
is ^ = co/L, \dipm\ < and / dz(^5'o -\- = 

to keep / dz^'^vE' = 1 + 0{S'i>^). Linearizing the cou- 
pled nonlinear-Schrodinger equations that describe the 



system we get /y 



ReG^, where x — ^lU: Gx 



y2((5^'i + (5^'_i)/VI, Gy ^ iV2((55'i - (5*_i)/\/I, and 



dt ^ 2dz^ ^ 



— gGy - ^(Gy + G*), 



where a — 2\ci\/L. To solve this equation we go to 
momentum space, a^{k) — J dzf^e'}cp{ikz) and b^{k) = 
J dzImG^ exp(iA:z), getting 
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Diagonalizing the matrix from Eq. ([6|) we see that there 
is instability for /a < 2 — q a.s in the Bogolubov spec- 
trum of this model [13]. Thus, the system is stable in 
the polar phase, and so small initial perturbations do 
not grow during the evolution towards broken-symmetry 
phase. The instability for g < 2 is responsible for the 
magnetization jump in Fig. [T]and the subsequent break- 
down of the linear approach (Fig. [Ud). 

To solve Eq. ^ with q{t) given by ^ we derive the 
equation for d^a^{t) / dt^ , keep leading order terms in the 
slow transition (tq ^ 1) and long-wavelength [k^ ja ^ 
2) fimits, and get that 



a-^{k,t) = Q;fcyAi(s)-|-/3fe;^Bi(s), 
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FIG. 4: (color) Dynamics of magnetic domains in fz. Black 
line (tq = 30), red line (rg = 50), green line (tq = 80). The 
arrow show direction of evolution on the main plot. Inset (a): 
early stages of £^{t) evolution. Inset (b): dependence of the 
typical post-transition domain size, ^, on quench time. The 
fit was done to tq > 30 data (see text for details). The figure 
shows results averaged over Nr — 44 runs; see [12 ] for units. 

where k = (q;^/2)^/'^, akx and Pkx are constants given 
by initial conditions, while Ai and Bi are Airy func- 
tions. From ([7]) we see that the instability arises from 
unbounded increase of the Bi(s) function happening for 
s > 0, i.e., fc^/a < 2 — g(t), which is a dynamical manifes- 
tation of the static result for unstable modes. This solu- 

"1/3 

tion works till t ^ t Tq when a significant increase of 
fx invalidates the linearized theory: this calculation rig- 
orously derives scaling ([3]) . Additionally, the solution ([7]) 
can be reliably used as long as tq 3> 1 or 37ms, which 
is also supported by numerics (Fig. [Iji). The quench 
time scale in the experiment [5| is much smaller than 
this bound. Finally, these results hold for any initial 
state spread over the k modes. 

The (re)scalings t/Tq^ and ^ ^ i ^ Tq'^ derived above 
in a ID system were also found by different means in a 2D 
spinor condensate |T^. A trivial extension of our mean- 
field analytical calculations to 2D and 3D systems shows 
that they hold for any number of spatial dimensions. 

To summarize, we have developed a theory of the dy- 
namics of symmetry-breaking in the quantum phase tran- 
sition inspired by the experiment [5] , but for the range of 
quench rates that are sufficiently slow so that the critical 
scalings can determine phase transition dynamics. This 
regime should be accessible by a "slower" version of the 
quench Our analysis points to a Kibble- Zurck-like 
scenario, where the state of the system departs from the 
old symmetric vacuum with a delay ~ t after the critical 
point was crossed. This sets up an initial post-transition 
state with a characteristic length scale ^ ([5]). This scale 
should determine the initial density of topological fea- 
tures. In our ID simulations textures appear, but we 
predict that in real 3D experiments other topological de- 
fects are created (as they were in Q), and the distance 
between them should be initially ^ ^. Such topological 
defects are more stable than textures so measurement of 



their density should be possible and would be a good test 
of the theory we have presented. 
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